################################################################################################
# This script takes the clean recoded SV data and uses the histogram from the Population Imm. 
# proposal to bootstrap the SV exp data. The reasoning behind this bootsrap is because it is  
# plausible that SV immigration proposals are a distinctly different sample
# from that observed in the QV immigration, Bootstrap will sample with replacement. 
# 
# It goes through 10000 bootstraps conditonal on imitating the Pop Imm Pref histogram and runs 
# SV Rules that would be used in the normal bootstrap for QV/SV. It saves the results of the
# rules and also a copy of the distribution preferences for each proposal in each bootstrap 
# iteration. Exports these files locally to the working directory.
# 
# This script takes the clean QV data and bootstraps  the exp data to determine welfare under 
# different rules to choose a class and cast the votes. The following rules are tested:
# 
# - Rule A: Choose class and cast vote as chosen in experiment
# - Rule B: 
# - Rule C: 
# - Rule D:  
# - Rule E:
# 
# Rewrites the Class and Vote columns based on these and calculates some welfare measures as in 
# Figures 4 and 5 and Table 2 of the Columbia Elections Paper.
# 
# All data is the saved in ./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/ 
# 
# Author: Luis S.
# Last Modified: 8.25.17
################################################################################################

library(dplyr)
library(gtools)

setwd("C:/Users/Luis/Dropbox/Research/SV vs QV California Experiment/Experiment Analysis & Results")

QV_data <- read.csv("California Stage 2 - SV & QV Data/Clean Data/California_MTurk_QV_Stage_2_clean_prefrecoded.csv", header = T, stringsAsFactors = F)
QV_data <- tibble::rownames_to_column(QV_data, var = "RowNum")
QV_data <- QV_data %>% mutate(RowNum = as.numeric(RowNum))

SV_data <- read.csv("California Stage 2 - SV & QV Data/Clean Data/California_MTurk_SV_Stage_2_clean_prefrecoded.csv", header = T, stringsAsFactors = F)

standardVoteWeights <- c(2,1.5,1.2,1)
# newVoteWeights <- standardVoteWeights + 1
# standardVoteWeights <- newVoteWeights

###################################
#### QV Vote & Welfare Summary ####
###################################

#--- Vote Summary and Welfare Summary Fns

QV_VoteSummary <- function(tempBootstrap)
{
  #--- Helper function for QV Summary
  
  detEffResult <- function(Position, Intensity)
  {
    # synthesize into a single dataframe and separate by positions
    temp <- data.frame(RefPosition = Position, PrefIntensity = Intensity, stringsAsFactors = F)
    
    # separate by infavor and against
    tempInFav <- temp %>% filter(RefPosition == "InFavor") 
    tempAgainst <- temp %>% filter(RefPosition == "Against")
    
    # compare intensity of preferences and determine winner
    tempInFav <- sum(tempInFav$PrefIntensity)
    tempAgainst <- sum(tempAgainst$PrefIntensity)
    effWinner <- ifelse(tempInFav == tempAgainst, sample(c("InFavor","Against"), 1, prob = c(.5,.5), replace = T),
                        ifelse(tempInFav > tempAgainst, "InFavor", "Against"))
    
    return(effWinner)
  }
  
  # #FOR TESTING ONLY
  # tempBootstrap <- BS_data
  # tempBootstrap <- expData
  
  #--- Replaces votes with actual weights 
  tempBootstrap[is.na(tempBootstrap)] <- 0
  
  for(row in 1:nrow(tempBootstrap))
  {
    # QV1
    if(tempBootstrap$QV1_Class[row] == "(1) Red") {voteWeight <- standardVoteWeights[1]}
    else if(tempBootstrap$QV1_Class[row] == "(2) Yellow") {voteWeight <- standardVoteWeights[2]} 
    else if(tempBootstrap$QV1_Class[row] == "(3) Green") {voteWeight <- standardVoteWeights[3]} 
    else {voteWeight <- standardVoteWeights[4]} 
    
    tempBootstrap$QV1_Vote_BilingualEduc[row] <- tempBootstrap$QV1_Vote_BilingualEduc[row] * voteWeight
    tempBootstrap$QV1_Vote_TeacherTenure[row] <- tempBootstrap$QV1_Vote_TeacherTenure[row] * voteWeight
    tempBootstrap$QV1_Vote_PublicBonds[row] <- tempBootstrap$QV1_Vote_PublicBonds[row] * voteWeight
    tempBootstrap$QV1_Vote_Immigration[row] <- tempBootstrap$QV1_Vote_Immigration[row] * voteWeight
    
    
    # QV2
    if(tempBootstrap$QV2_Class[row] == "(1) Red") {voteWeight <- standardVoteWeights[1]}
    else if(tempBootstrap$QV2_Class[row] == "(2) Yellow") {voteWeight <- standardVoteWeights[2]} 
    else if(tempBootstrap$QV2_Class[row] == "(3) Green") {voteWeight <- standardVoteWeights[3]} 
    else {voteWeight <- standardVoteWeights[4]} 
    
    tempBootstrap$QV2_Vote_BilingualEduc[row] <- tempBootstrap$QV2_Vote_BilingualEduc[row] * voteWeight
    tempBootstrap$QV2_Vote_TeacherTenure[row] <- tempBootstrap$QV2_Vote_TeacherTenure[row] * voteWeight
    tempBootstrap$QV2_Vote_PublicBonds[row] <- tempBootstrap$QV2_Vote_PublicBonds[row] * voteWeight
    tempBootstrap$QV2_Vote_Immigration[row] <- tempBootstrap$QV2_Vote_Immigration[row] * voteWeight
  }
  
  # Determine count w/o QV
  tempRef1 <- group_by(tempBootstrap, Ref1_BilingualEduc) %>% summarise(BilingualEduc = n())
  tempRef2 <- group_by(tempBootstrap, Ref2_TeacherTenure) %>% summarise(TeacherTenure = n())
  tempRef3 <- group_by(tempBootstrap, Ref3_PublicBonds) %>% summarise(PublicBonds = n())
  tempRef4 <- group_by(tempBootstrap, Ref4_Immigration) %>% summarise(Immigration = n())
  
  tempQV_noQV <- bind_cols(tempRef1,tempRef2,tempRef3,tempRef4) %>% select(-matches("Ref"))
  
  # Determine Weighted vote count & add vote side w/o QV  
  
  tempRef1 <- group_by(tempBootstrap, Ref1_BilingualEduc) %>% summarise(BilingualEduc_QV1 = sum(QV1_Vote_BilingualEduc), BilingualEduc_QV2 = sum(QV2_Vote_BilingualEduc))
  tempRef2 <- group_by(tempBootstrap, Ref2_TeacherTenure) %>% summarise(TeacherTenure_QV1 = sum(QV1_Vote_TeacherTenure), TeacherTenure_QV2 = sum(QV2_Vote_TeacherTenure))
  tempRef3 <- group_by(tempBootstrap, Ref3_PublicBonds) %>% summarise(PublicBonds_QV1 = sum(QV1_Vote_PublicBonds), PublicBonds_QV2 = sum(QV2_Vote_PublicBonds))
  tempRef4 <- group_by(tempBootstrap, Ref4_Immigration) %>% summarise(Immigration_QV1 = sum(QV1_Vote_Immigration), Immigration_QV2 = sum(QV2_Vote_Immigration))
  
  QV_summary <- bind_cols(tempQV_noQV,tempRef1,tempRef2,tempRef3,tempRef4) %>% select(-matches("Ref"))
  
  # #### Add Simple Votes to QV votes
  # QV_summary_old <- QV_summary
  # QV_summary <- QV_summary %>% mutate(BilingualEduc_QV1 = BilingualEduc_QV1 + BilingualEduc,
  #                                     BilingualEduc_QV2 = BilingualEduc_QV2 + BilingualEduc,
  #                                     TeacherTenure_QV1 = TeacherTenure_QV1 + TeacherTenure,
  #                                     TeacherTenure_QV2 = TeacherTenure_QV2 + TeacherTenure,
  #                                     PublicBonds_QV1 = PublicBonds_QV1 + PublicBonds,
  #                                     PublicBonds_QV2 = PublicBonds_QV2 + PublicBonds,
  #                                     Immigration_QV1 = Immigration_QV1 + Immigration,
  #                                     Immigration_QV2 = Immigration_QV2 + Immigration)
  # 
  # #### End
  
  
  QV_summary <- as.data.frame(t(QV_summary))
  
  QV_summary <- tibble::rownames_to_column(QV_summary, "Position") %>% arrange(Position)
  
  names(QV_summary) <- c("Position", "Abstain", "Against", "InFavor")
  row.names(QV_summary) <- seq(1,nrow(QV_summary),1)
  
  
  assign("Testing_QV_Summary_Pre", QV_summary, envir=.GlobalEnv)
  
  # Calculate desired values for each of the referenda 
  
  for(ref in 1:4) # iterates in multiples of four to get data for each referenda
  {
    index_noQV <- (ref*3)-2
    index_QV1 <- (ref*3)-1
    index_QV2 <- (ref*3)
    
    # Calculating simple majority result w/o QV (rows 1,4,7,10) & Determine the simple majority winner
    if(QV_summary$InFavor[index_noQV] == QV_summary$Against[index_noQV]) 
    {
      QV_summary$NoQV_Majority[index_noQV] <- "Tie"
      QV_summary$NoQV_Majority[index_QV1] <- "Tie"
      QV_summary$NoQV_Majority[index_QV2] <- "Tie"
      
      tieBreaker <- ifelse(sample(2, 1, replace=T) > 1, "InFavor", "Against")
      
      QV_summary$NoQV_Result[index_noQV] <- tieBreaker
      QV_summary$NoQV_Result[index_QV1] <- tieBreaker
      QV_summary$NoQV_Result[index_QV2] <- tieBreaker
      
    } else if(QV_summary$InFavor[index_noQV] > QV_summary$Against[index_noQV]) 
    {
      QV_summary$NoQV_Majority[index_noQV] <- "InFavor"
      QV_summary$NoQV_Majority[index_QV1] <- "InFavor"
      QV_summary$NoQV_Majority[index_QV2] <- "InFavor"
    } else {
      QV_summary$NoQV_Majority[index_noQV] <- "Against"
      QV_summary$NoQV_Majority[index_QV1] <- "Against"
      QV_summary$NoQV_Majority[index_QV2] <- "Against"
    }
    
    QV_summary <- mutate(QV_summary, NoQV_Result = ifelse(NoQV_Majority == "Tie", NoQV_Result, NoQV_Majority))
    
    
    # Calculating QV1 result (rows 2,5,8,11) & Determine the QV1 winner
    if(QV_summary$InFavor[index_QV1] == QV_summary$Against[index_QV1]) 
    {
      QV_summary$QV1_Majority[index_noQV] <- "Tie"
      QV_summary$QV1_Majority[index_QV1] <- "Tie"
      QV_summary$QV1_Majority[index_QV2] <- "Tie"
      
      tieBreaker <- ifelse(sample(2, 1, replace=T) > 1, "InFavor", "Against")
      
      QV_summary$QV1_Result[index_noQV] <- tieBreaker
      QV_summary$QV1_Result[index_QV1] <- tieBreaker
      QV_summary$QV1_Result[index_QV2] <- tieBreaker
      
    } else if(QV_summary$InFavor[index_QV1] > QV_summary$Against[index_QV1]) 
    {
      QV_summary$QV1_Majority[index_noQV] <- "InFavor"
      QV_summary$QV1_Majority[index_QV1] <- "InFavor"
      QV_summary$QV1_Majority[index_QV2] <- "InFavor"
    } else {
      QV_summary$QV1_Majority[index_noQV] <- "Against"
      QV_summary$QV1_Majority[index_QV1] <- "Against"
      QV_summary$QV1_Majority[index_QV2] <- "Against"
    }
    
    QV_summary <- mutate(QV_summary, QV1_Result = ifelse(QV1_Majority == "Tie", QV1_Result, QV1_Majority))
    
    # Calculating QV2 result (rows 2,5,8,11) & Determine the QV2 winner
    if(QV_summary$InFavor[index_QV2] == QV_summary$Against[index_QV2]) 
    {
      QV_summary$QV2_Majority[index_noQV] <- "Tie"
      QV_summary$QV2_Majority[index_QV1] <- "Tie"
      QV_summary$QV2_Majority[index_QV2] <- "Tie"
      
      tieBreaker <- ifelse(sample(2, 1, replace=T) > 1, "InFavor", "Against")
      
      QV_summary$QV2_Result[index_noQV] <- tieBreaker
      QV_summary$QV2_Result[index_QV1] <- tieBreaker
      QV_summary$QV2_Result[index_QV2] <- tieBreaker
      
    } else if(QV_summary$InFavor[index_QV2] > QV_summary$Against[index_QV2]) 
    {
      QV_summary$QV2_Majority[index_noQV] <- "InFavor"
      QV_summary$QV2_Majority[index_QV1] <- "InFavor"
      QV_summary$QV2_Majority[index_QV2] <- "InFavor"
    } else {
      QV_summary$QV2_Majority[index_noQV] <- "Against"
      QV_summary$QV2_Majority[index_QV1] <- "Against"
      QV_summary$QV2_Majority[index_QV2] <- "Against"
    }
    
    QV_summary <- mutate(QV_summary, QV2_Result = ifelse(QV2_Majority == "Tie", QV2_Result, QV2_Majority))
  }
  
  # clean summary table
  QV_summary <- QV_summary[grep("[^2])*$", QV_summary$Position), ]
  QV_summary <- QV_summary[grep("[^1])*$", QV_summary$Position), ]
  
  assign("Testing_QV_Summary_Mid", QV_summary, envir=.GlobalEnv)
  
  QV_summary <- select(QV_summary, Position, matches("Result")) %>% rename(Referenda = Position) %>% mutate(Eff_Result = NA)
  
  # Determine the efficient outcome and add it to QV Summary
  
  QV_summary$Eff_Result[QV_summary$Referenda == "BilingualEduc"] <- detEffResult(tempBootstrap$Ref1_BilingualEduc, tempBootstrap$EducationPref)
  QV_summary$Eff_Result[QV_summary$Referenda == "Immigration"] <- detEffResult(tempBootstrap$Ref4_Immigration, tempBootstrap$ImmigrationPref)
  QV_summary$Eff_Result[QV_summary$Referenda == "PublicBonds"] <- detEffResult(tempBootstrap$Ref3_PublicBonds, tempBootstrap$BondsPref)
  QV_summary$Eff_Result[QV_summary$Referenda == "TeacherTenure"] <- detEffResult(tempBootstrap$Ref2_TeacherTenure, tempBootstrap$TeachersPref)
  
  # Determine if the BV helped the minority win
  
  QV_summary <- QV_summary %>%
    mutate(QV1_Winner = ifelse(NoQV_Result == QV1_Result, "Majority", "Minority"),
           QV2_Winner = ifelse(NoQV_Result == QV2_Result, "Majority", "Minority")) %>%
    mutate(IsEff_NQV = ifelse(NoQV_Result == Eff_Result, "Efficient", "Inefficient"),
           IsEff_QV1 = ifelse(QV1_Result == Eff_Result, "Efficient", "Inefficient"),
           IsEff_QV2 = ifelse(QV2_Result == Eff_Result, "Efficient", "Inefficient"))
  
  assign("QV_summary", QV_summary, envir=.GlobalEnv)
  
  NQV_Results <- QV_summary$NoQV_Result
  Eff_Results <- QV_summary$Eff_Result
  QV1_Results <- QV_summary$QV1_Result
  QV2_Results <- QV_summary$QV2_Result
  QV1_Winner <- QV_summary$QV1_Winner
  QV2_Winner <- QV_summary$QV2_Winner
  
  return(c(NQV_Results,Eff_Results,QV1_Results,QV2_Results,QV1_Winner,QV2_Winner))
}


QV_WelfareSummary <- function(tempBootstrap, majorityResults)
{
  # We want to calculate the aggregate welfare of the QV mechanism vs the simple majority
  # Thus, we want to see what the total welfare over the four elections is if we allow
  # the majority to be as in simple majority, versus as if we were in the QV world. Note 
  # that by majority we mean the winning side
  
  # QV_summary contains the side of the majority with and without QV votes. Rows 1,2,3,4 
  # correspond to BilingualEduc, Immigration, PublicBonds, and TeachersTenure, respectively. 
  
  # Aggregate Welfare under Simple Majority
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == QV_summary$NoQV_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == QV_summary$NoQV_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == QV_summary$NoQV_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == QV_summary$NoQV_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_noQV <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  # Aggregate Welfare under Efficiency
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == QV_summary$Eff_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == QV_summary$Eff_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == QV_summary$Eff_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == QV_summary$Eff_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_Eff <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  
  # Aggregate Welfare under Storable Votes - QV1
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == QV_summary$QV1_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == QV_summary$QV1_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == QV_summary$QV1_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == QV_summary$QV1_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_QV1 <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  # Aggregate Welfare under Storable Votes - QV2
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == QV_summary$QV2_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == QV_summary$QV2_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == QV_summary$QV2_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == QV_summary$QV2_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_QV2 <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  return(c(totalWelfare_noQV, totalWelfare_Eff, totalWelfare_QV1, totalWelfare_QV2))
}

#--- Function to store Vote and Welfare Summaries

saveSummary <- function(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
{
  # Saves what side won in simple majority
  WelfareResults$NQV_Result_BilingualEduc[bootstrapNum] <- voteSummary[1]
  WelfareResults$NQV_Result_Immigration[bootstrapNum] <- voteSummary[2]
  WelfareResults$NQV_Result_PublicBonds[bootstrapNum] <- voteSummary[3]
  WelfareResults$NQV_Result_TeacherTenure[bootstrapNum] <- voteSummary[4]
  
  # Saves what side should have won under efficient outcomes
  WelfareResults$Eff_Result_BilingualEduc[bootstrapNum] <- voteSummary[5]
  WelfareResults$Eff_Result_Immigration[bootstrapNum] <- voteSummary[6]
  WelfareResults$Eff_Result_PublicBonds[bootstrapNum] <- voteSummary[7]
  WelfareResults$Eff_Result_TeacherTenure[bootstrapNum] <- voteSummary[8]
  
  # Save what side won under QV
  WelfareResults$QV1_Result_BilingualEduc[bootstrapNum] <- voteSummary[9]
  WelfareResults$QV1_Result_Immigration[bootstrapNum] <- voteSummary[10]
  WelfareResults$QV1_Result_PublicBonds[bootstrapNum] <- voteSummary[11]
  WelfareResults$QV1_Result_TeacherTenure[bootstrapNum] <- voteSummary[12]
  
  WelfareResults$QV2_Result_BilingualEduc[bootstrapNum] <- voteSummary[13]
  WelfareResults$QV2_Result_Immigration[bootstrapNum] <- voteSummary[14]
  WelfareResults$QV2_Result_PublicBonds[bootstrapNum] <- voteSummary[15]
  WelfareResults$QV2_Result_TeacherTenure[bootstrapNum] <- voteSummary[16]
  
  # Save whether the minority or majority won (assumes that majority is simple majority)
  WelfareResults$QV1_Winner_BilingualEduc[bootstrapNum] <- voteSummary[17]
  WelfareResults$QV1_Winner_Immigration[bootstrapNum] <- voteSummary[18]
  WelfareResults$QV1_Winner_PublicBonds[bootstrapNum] <- voteSummary[19]
  WelfareResults$QV1_Winner_TeacherTenure[bootstrapNum] <- voteSummary[20]
  
  WelfareResults$QV2_Winner_BilingualEduc[bootstrapNum] <- voteSummary[21]
  WelfareResults$QV2_Winner_Immigration[bootstrapNum] <- voteSummary[22]
  WelfareResults$QV2_Winner_PublicBonds[bootstrapNum] <- voteSummary[23]
  WelfareResults$QV2_Winner_TeacherTenure[bootstrapNum] <- voteSummary[24]
  
  # save welfare summary 
  
  WelfareResults$AggregateWelfare_noQV[bootstrapNum] <- welfareSummary[1]
  WelfareResults$AggregateWelfare_Eff[bootstrapNum] <- welfareSummary[2]
  WelfareResults$AggregateWelfare_QV1[bootstrapNum] <- welfareSummary[3]
  WelfareResults$AggregateWelfare_QV2[bootstrapNum] <- welfareSummary[4]
  
  return(WelfareResults)
}


################################
#### Vote Casting Functions ####
###############################%

#--- helper functions

refsNotAbstained <- function(subjectData)
{
  refs <- c()
  
  if(subjectData$Ref1_BilingualEduc != "Abstain"){ refs <- c(refs, "BilingualEduc")}
  if(subjectData$Ref2_TeacherTenure != "Abstain"){ refs <- c(refs, "TeacherTenure")}
  if(subjectData$Ref3_PublicBonds != "Abstain"){ refs <- c(refs, "PublicBonds")}
  if(subjectData$Ref4_Immigration != "Abstain"){ refs <- c(refs, "Immigration")}
  
  return(refs) 
}

detVoteVector <- function(chosenRefs)
{
  voteVector <- data.frame(BilingualEduc = NA, TeacherTenure = NA, PublicBonds = NA, Immigration = NA)
  for(issue in chosenRefs) { voteVector[[issue]][1] <- 1}
  return(voteVector)
}

#--- casting functions

castStatModel <- function(subjectData, validRefs, rho, gamma, eta)
{
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  
  #--- choose a valid vote class 
  subjectPrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  subjectPrefStrength <- sort(subjectPrefStrength, decreasing = T)
  
  chosenClass <- ifelse(subjectPrefStrength[1] > rho*subjectPrefStrength[2], "(1) Red",
                        ifelse(subjectPrefStrength[2] > gamma*subjectPrefStrength[3], "(2) Yellow",
                               ifelse(subjectPrefStrength[3] > rho*subjectPrefStrength[4], "(3) Green", "(4) Blue")))
  
  #--- cast on the issues with most points (excl. abstain)
  numChosenRefs <- as.numeric(regmatches(chosenClass, regexpr("\\d+", chosenClass)))
  
  # match issues in which subject did not abstain with corresponding pref strength
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  
  validRefIndex <- viableReferenda %in% validRefs
  viableReferenda <- viableReferenda[validRefIndex]
  viablePrefStrength <- viablePrefStrength[validRefIndex]
  
  chosenRefs <- c()
  
  # find index of next max and choose it.
  for(iter in 1:numChosenRefs) 
  {
    nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
    nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
    
    chosenRefs <- c(chosenRefs, viableReferenda[nextRefIndex])
    
    validRefIndex <- !viableReferenda %in% chosenRefs
    viableReferenda <- viableReferenda[validRefIndex]
    viablePrefStrength <- viablePrefStrength[validRefIndex]
  }
  
  results <- list()
  results[[1]] <- chosenClass 
  results[[2]] <- chosenRefs
  
  return(results)
}

castStatModel_classAndVote <- function(subjectData, rho, gamma, eta, eps, mu)
{
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  # rho = 1.52; gamma <- 1.19; eta <- 1.15; eps <- .51; mu <- .2 
  # 
  
  # set.seed(7)
  
  #--- choose a valid vote class or make a mistake
  
  classIsCorrect <- ifelse(runif(1,0,1) <= eps, F, T) # determine if subject makes a mistake or not
  
  
  # find correct class
  
  subjectPrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  subjectPrefStrength <- sort(subjectPrefStrength, decreasing = T)
  
  
  tieBreaker <- round(runif(1,0,1)) + 1
  indiff_RedYellow <- c("(1) Red", "(2) Yellow")
  indiff_YellowGreen <- c("(2) Yellow", "(3) Green")
  indiff_GreenBlue <- c("(3) Green", "(4) Blue")
  
  correctClass <- ifelse(subjectPrefStrength[1] > rho*subjectPrefStrength[2], "(1) Red",
                         ifelse(subjectPrefStrength[1] == rho*subjectPrefStrength[2], indiff_RedYellow[tieBreaker],
                                ifelse(subjectPrefStrength[2] > gamma*subjectPrefStrength[3], "(2) Yellow",
                                       ifelse(subjectPrefStrength[2] == gamma*subjectPrefStrength[3], indiff_YellowGreen[tieBreaker],
                                              ifelse(subjectPrefStrength[3] > eta*subjectPrefStrength[4], "(3) Green",
                                                     ifelse(subjectPrefStrength[3] == eta*subjectPrefStrength[4], indiff_GreenBlue[tieBreaker], "(4) Blue"))))))
  
  
  # save either correct class or mistaken class
  
  if(classIsCorrect) { chosenClass <- correctClass } else {
    
    # make sure not to choose class that will force vote on abstain       ### Not true... 
    voteClasses <- c("(1) Red","(2) Yellow", "(3) Green", "(4) Blue")
    # voteClasses <- voteClasses[1:numValidRefs]
    
    # eliminate the correct vote class
    validClasses <- voteClasses[!voteClasses %in% correctClass]
    chosenClass <- sample(x = validClasses, size = 1, prob = rep(1,length(validClasses)), rep = T)
  } 
  
  
  #--- cast votes either monotonically or non-monotonically
  
  numChosenRefs <- as.numeric(regmatches(chosenClass, regexpr("\\d+", chosenClass)))
  
  # determine if vote is monotone by default (blue class) else toss coin to decide
  
  if(numChosenRefs == 4)
  {
    chosenRefs <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  } else {
    
    voteIsMonotone <- ifelse(runif(1,0,1) <= mu, F, T) # determine if subject makes a mistake or not
    
    if(voteIsMonotone)
    {
      currentChoice <- areMonotone[[subjectData$RowNum]][[numChosenRefs]] 
      if(nrow(currentChoice) == 0) {currentChoice <- areNotMonotone[[subjectData$RowNum]][[numChosenRefs]] } # if the choice is empty then go for the nonmonotone
      chosenRefs <- unlist(currentChoice[sample(seq(1:nrow(currentChoice)),1),]) 
    } else {
      currentChoice <- areNotMonotone[[subjectData$RowNum]][[numChosenRefs]] 
      if(nrow(currentChoice) == 0) {currentChoice <- areMonotone[[subjectData$RowNum]][[numChosenRefs]] } # if the choice is empty then go for monotone
      chosenRefs <- unlist(currentChoice[sample(seq(1:nrow(currentChoice)),1),]) 
    }
  }
  
  results <- list()
  results[[1]] <- chosenClass 
  results[[2]] <- chosenRefs
  
  return(results)
}


castNormativeModel <- function(subjectData, validRefs)
{
  # row <- 2
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  
  voteClasses <- c("(1) Red","(2) Yellow", "(3) Green", "(4) Blue")
  
  #--- choose a valid vote class
  subjectPrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  
  proportionalityFactor <- 2/sqrt(sum(subjectPrefStrength^2))
  
  normativeVoteWeights <- subjectPrefStrength*proportionalityFactor
  normativeVoteWeights <- sort(normativeVoteWeights, decreasing = T)
  
  # minimize choose the vector with the least euclidean distance
  calcEuclidDist <- function(v1,v2) {return(sqrt(sum((v2-v1)^2)))}
  
  redVector <- c(2,0,0,0)
  yellowVector <- c(1.5,1.5,0,0)
  greenVector <- c(1.2,1.2,1.2,0)
  blueVector <- c(1,1,1,1)
  
  subjectDistanceScores <- c(calcEuclidDist(normativeVoteWeights, redVector),
                             calcEuclidDist(normativeVoteWeights, yellowVector),
                             calcEuclidDist(normativeVoteWeights, greenVector),
                             calcEuclidDist(normativeVoteWeights, blueVector))
  
  subjectDistanceScores <- subjectDistanceScores[1:length(validRefs)] # select only the vote classes which will not force vote on abstain
  
  chosenClass <- which(min(subjectDistanceScores) == subjectDistanceScores) # choose the one that minimizes the distance
  chosenClass <- voteClasses[chosenClass]
  
  #--- cast on the issues with most points (excl. abstain)
  numChosenRefs <- as.numeric(regmatches(chosenClass, regexpr("\\d+", chosenClass)))
  
  # match issues in which subject did not abstain with corresponding pref strength
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  
  validRefIndex <- viableReferenda %in% validRefs
  viableReferenda <- viableReferenda[validRefIndex]
  viablePrefStrength <- viablePrefStrength[validRefIndex]
  
  chosenRefs <- c()
  
  # find index of next max and choose it.
  for(iter in 1:numChosenRefs) 
  {
    nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
    nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
    
    chosenRefs <- c(chosenRefs, viableReferenda[nextRefIndex])
    
    validRefIndex <- !viableReferenda %in% chosenRefs
    viableReferenda <- viableReferenda[validRefIndex]
    viablePrefStrength <- viablePrefStrength[validRefIndex]
  }
  
  results <- list()
  results[[1]] <- chosenClass 
  results[[2]] <- chosenRefs
  
  return(results)
}

castRand <- function(subjectData, validRefs, ruleFreq = c(1,1,1,1), castOnHighest = F)
{
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  numValidRefs <- length(validRefs)
  voteClasses <- c("(1) Red","(2) Yellow", "(3) Green", "(4) Blue")
  ruleFreq <- ruleFreq[1:numValidRefs]
  
  # choose a valid vote class
  chosenClass <- sample(1:numValidRefs, 1, prob = ruleFreq, replace = T)
  chosenClass <- voteClasses[chosenClass]
  
  # cast votes
  numChosenRefs <- as.numeric(regmatches(chosenClass, regexpr("\\d+", chosenClass)))
  
  if(castOnHighest) # votes are cast on the issues with most points (excl. abstain)
  {
    # match issues in which subject did not abstain with corresponding pref strength
    viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
    viablePrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
    
    validRefIndex <- viableReferenda %in% validRefs
    viableReferenda <- viableReferenda[validRefIndex]
    viablePrefStrength <- viablePrefStrength[validRefIndex]
    
    chosenRefs <- c()
    
    # find index of next max and choose it.
    for(iter in 1:numChosenRefs) 
    {
      nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
      nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
      
      chosenRefs <- c(chosenRefs, viableReferenda[nextRefIndex])
      
      validRefIndex <- !viableReferenda %in% chosenRefs
      viableReferenda <- viableReferenda[validRefIndex]
      viablePrefStrength <- viablePrefStrength[validRefIndex]
    }
  } else { chosenRefs <- sample(validRefs, numChosenRefs, prob = rep(1,numValidRefs), replace = F) } # votes are cast randomly (excl. abstain)
  
  results <- list()
  results[[1]] <- chosenClass 
  results[[2]] <- chosenRefs
  
  return(results)
}


chooseNormativeClass <- function(subjectData, validRefs)
{
  # row <- 2
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  
  voteClasses <- c("(1) Red","(2) Yellow", "(3) Green", "(4) Blue")
  
  #--- choose a valid vote class
  subjectPrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  
  proportionalityFactor <- 2/sqrt(sum(subjectPrefStrength^2))
  
  normativeVoteWeights <- subjectPrefStrength*proportionalityFactor
  normativeVoteWeights <- sort(normativeVoteWeights, decreasing = T)
  
  # minimize choose the vector with the least euclidean distance
  calcEuclidDist <- function(v1,v2) {return(sqrt(sum((v2-v1)^2)))}
  
  redVector <- c(standardVoteWeights[1],0,0,0)
  yellowVector <- c(standardVoteWeights[2],standardVoteWeights[2],0,0)
  greenVector <- c(standardVoteWeights[3],standardVoteWeights[3],standardVoteWeights[3],0)
  blueVector <- rep(standardVoteWeights[4],4)
  
  subjectDistanceScores <- c(calcEuclidDist(normativeVoteWeights, redVector),
                             calcEuclidDist(normativeVoteWeights, yellowVector),
                             calcEuclidDist(normativeVoteWeights, greenVector),
                             calcEuclidDist(normativeVoteWeights, blueVector))
  
  subjectDistanceScores <- subjectDistanceScores[1:length(validRefs)] # select only the vote classes which will not force vote on abstain
  
  chosenClass <- which(min(subjectDistanceScores) == subjectDistanceScores) # choose the one that minimizes the distance
  chosenClass <- voteClasses[chosenClass]
  
  return(chosenClass)
  
}

chooseRandClass <- function(subjectData, validRefs)
{
  ruleFreq = c(1,1,1,1)
  
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  numValidRefs <- length(validRefs)
  voteClasses <- c("(1) Red","(2) Yellow", "(3) Green", "(4) Blue")
  ruleFreq <- ruleFreq[1:numValidRefs]
  
  # choose a valid vote class
  chosenClass <- sample(1:numValidRefs, 1, prob = ruleFreq, replace = T)
  chosenClass <- voteClasses[chosenClass]
  
  return(chosenClass)
}


castVotesMonotically <- function(subjectData, validRefs, chosenClass)
{
  # row <- 2
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  
  voteClasses <- c("(1) Red","(2) Yellow", "(3) Green", "(4) Blue")
  
  #--- cast on the issues with most points (excl. abstain)
  numChosenRefs <- as.numeric(regmatches(chosenClass, regexpr("\\d+", chosenClass)))
  
  # match issues in which subject did not abstain with corresponding pref strength
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(subjectData$EducationPref[1], subjectData$TeachersPref[1], subjectData$BondsPref[1], subjectData$ImmigrationPref[1])
  
  validRefIndex <- viableReferenda %in% validRefs
  viableReferenda <- viableReferenda[validRefIndex]
  viablePrefStrength <- viablePrefStrength[validRefIndex]
  
  chosenRefs <- c()
  
  # find index of next max and choose it.
  for(iter in 1:numChosenRefs) 
  {
    nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
    nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
    
    chosenRefs <- c(chosenRefs, viableReferenda[nextRefIndex])
    
    validRefIndex <- !viableReferenda %in% chosenRefs
    viableReferenda <- viableReferenda[validRefIndex]
    viablePrefStrength <- viablePrefStrength[validRefIndex]
  }
  
  results <- list()
  results[[1]] <- chosenClass 
  results[[2]] <- chosenRefs
  
  return(results)
}


castVotesRandomly <- function(subjectData, validRefs, chosenClass)
{
  # subjectData <- expData[row,]
  # validRefs <- refsNotAbstained(expData[row,])
  
  numValidRefs <- length(validRefs)
  
  # cast votes
  numChosenRefs <- as.numeric(regmatches(chosenClass, regexpr("\\d+", chosenClass)))
  
  chosenRefs <- sample(validRefs, numChosenRefs, prob = rep(1,numValidRefs), replace = F)
  
  results <- list()
  results[[1]] <- chosenClass 
  results[[2]] <- chosenRefs
  
  return(results)
}



###################################
##### QV Rule A - AS IS Model #####
##################################%

ruleA <- function(expData, WelfareResults, bootstrapNum)
{
  # This rule uses chooses the voting class and cast the votes as they were cast 
  # by the original subjects. Thus no modification is needed and we just run the welfare measures.
  
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleAResults
  # bootstrapNum <- 1
  
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleAResults", WelfareResults, envir=.GlobalEnv)
}


##################################
##### QV Rule B - Stat Model #####
#################################%

# Use estimated rho, gamma, and eta (corresponding to the pairwise ratios of the issues ordered by pref strength)
# which minimize the number of misclassified individuals (according to vote classes) to choose the vote classes.
# Then, cast vote on most preferred issue

StatModelParams <- c(1.25,1.05,1)

ruleB <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleAResults
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    # QV1
    results <- castStatModel(expData[row,], validRefs, StatModelParams[1], StatModelParams[2], StatModelParams[3])
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castStatModel(expData[row,], validRefs, StatModelParams[1], StatModelParams[2], StatModelParams[3])
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleBResults", WelfareResults, envir=.GlobalEnv)
}



##########################################
##### QV Rule B - Model 2 with noise #####
#########################################%

# Use estimated rho, gamma, and eta (corresponding to the pairwise ratios of the issues ordered by pref strength)
# along with epsilon and mu which correspond to the prob of making a mistake in choosing a class and casting monotonically respectively.  

# rho, gamma, eta, eps, mu

estimatedParams_QV1 <- c(1.329,1.186,1.388,.498,.211)
estimatedParams_QV2 <- c(1.251,1.180,1.452,.492,.262)


#--- Load estimated 

QV_AllPos <- read.csv("California Stage 2 - SV & QV Data/Clean Data/California_MTurk_QV_Stage_2_clean&Recoded_AllPos.csv", header = T, stringsAsFactors = F)
QV_AllPos[is.na(QV_AllPos)] <- 0

IsVoteMonotone <- function(candidateRefs, Prefs, Refs)
{
  
  indexes <- c()
  
  # find the pref intensities that correspond to candidate refs 
  for(i in 1:length(candidateRefs)) { indexes <- c(indexes,which(Refs == candidateRefs[i])) }
  candidatePrefs <- Prefs[indexes]
  
  # sort candidate prefs and refs in decreasing order
  candidateRefs <- candidateRefs[order(candidatePrefs, decreasing = T)]
  candidatePrefs <- candidatePrefs[order(candidatePrefs, decreasing = T)]
  
  # iterate through the candidates and check that there is no issue with more points that wasnt chosen
  maxChecks <- length(candidateRefs)
  isMonotone <- T
  
  for(iter in 1:maxChecks)
  {
    currentPref <- candidatePrefs[iter] 
    currentRef <- candidateRefs[iter]
    
    if(max(Prefs) > currentPref) { isMonotone <- F } else {
      # remove current ones from consideration
      
      alreadyCheckedIndex <- which(Refs == currentRef)
      Prefs <- Prefs[-alreadyCheckedIndex]
      Refs <- Refs[-alreadyCheckedIndex]
    }
    
    if(!isMonotone) { break()} # once you find a monotonicity mistake, isMonotone == F, exit loop
    
  }
  
  return(isMonotone)
}

# find all posible ways to cast the votes for a given rule
RefList <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration") 

feasibleRed <- combinations(r = 1, n = 4, v = RefList)
feasibleYellow <- combinations(r = 2, n = 4, v = RefList)
feasibleGreen <- combinations(r = 3, n = 4, v = RefList)

areMonotone <- list()
areNotMonotone <- list()

for(row in 1:nrow(QV_AllPos))
{
  
  #--- Determine number of ways in which subject could make monotonicity mistake under feasible classes
  
  # find number of abstains in case we want to restrict epsilon to only those that dont force a vote on abstain
  voteDirection <- c(QV_AllPos$Ref1_BilingualEduc[row], QV_AllPos$Ref2_TeacherTenure[row], QV_AllPos$Ref3_PublicBonds[row], QV_AllPos$Ref4_Immigration[row]) #must match ref list order
  RefList <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration") # must match vote direction order
  PrefStrList <- c(QV_AllPos$EducationPref[row], QV_AllPos$TeachersPref[row], QV_AllPos$BondsPref[row], QV_AllPos$ImmigrationPref[row]) # must match order of two above
  
  
  # create vectors to store results
  monotoneRed <- data.frame()
  monotoneYellow <- data.frame()
  monotoneGreen <- data.frame()
  
  notMonotoneRed <- data.frame()
  notMonotoneYellow <- data.frame()
  notMonotoneGreen <- data.frame()
  
  # iterate and check
  for(candidate in 1:nrow(feasibleRed)) 
  { 
    if(IsVoteMonotone(feasibleRed[candidate,], PrefStrList, RefList)) { monotoneRed <- bind_rows(monotoneRed, as.data.frame(t(feasibleRed[candidate,]),stringsAsFactors = F)) } else { notMonotoneRed <- bind_rows(notMonotoneRed, as.data.frame(t(feasibleRed[candidate,]),stringsAsFactors = F)) }
  }
  
  for(candidate in 1:nrow(feasibleYellow)) 
  { 
    if(IsVoteMonotone(feasibleYellow[candidate,], PrefStrList, RefList)) { monotoneYellow <- bind_rows(monotoneYellow, as.data.frame(t(feasibleYellow[candidate,]),stringsAsFactors = F) ) } else { notMonotoneYellow <- bind_rows(notMonotoneYellow, as.data.frame(t(feasibleYellow[candidate,]),stringsAsFactors = F)) }  
  }
  
  for(candidate in 1:nrow(feasibleGreen)) 
  { 
    if(IsVoteMonotone(feasibleGreen[candidate,], PrefStrList, RefList)) { monotoneGreen <- bind_rows(monotoneGreen, as.data.frame(t(feasibleGreen[candidate,]),stringsAsFactors = F) ) } else { notMonotoneGreen <- bind_rows(notMonotoneGreen, as.data.frame(t(feasibleGreen[candidate,]),stringsAsFactors = F)) }
  }
  
  
  #--- Save the output in lists
  
  areMonotone[[row]] <- list(monotoneRed, monotoneYellow, monotoneGreen, QV_AllPos$MID[row])
  areNotMonotone[[row]] <- list(notMonotoneRed, notMonotoneYellow, notMonotoneGreen, QV_AllPos$MID[row])
  
}



#--- Estimate Model

ruleB_classAndVote <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleBResults_classAndVote
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    
    # QV1
    results <- castStatModel_classAndVote(expData[row,], estimatedParams_QV1[1], estimatedParams_QV1[2], estimatedParams_QV1[3], estimatedParams_QV1[4], estimatedParams_QV1[5])
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castStatModel_classAndVote(expData[row,], estimatedParams_QV2[1], estimatedParams_QV2[2], estimatedParams_QV2[3], estimatedParams_QV2[4], estimatedParams_QV2[5])
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  
  # get summaries and store them
  
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleBResults_classAndVote", WelfareResults, envir=.GlobalEnv)
}



#######################################
##### QV Rule C - Normative Model #####
######################################%

ruleC <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleAResults
  # bootstrapNum <- 1
  
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    # QV1
    results <- castNormativeModel(expData[row,], validRefs)
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castNormativeModel(expData[row,], validRefs)
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleCResults", WelfareResults, envir=.GlobalEnv)
}





###############################################################
#### QV Rule D (Uniform Vote Class) - Cast Vote on Highest ####
##############################################################%

# Want to choose a random class (using a uniform distribution over vote classes) and then allocate votes on highest ranked refs.

ruleD_UniformVoteClass <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleDResults
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    # QV1
    results <- castRand(expData[row,], validRefs, ruleFreq = c(1,1,1,1), castOnHighest = T)
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castRand(expData[row,], validRefs, ruleFreq = c(1,1,1,1), castOnHighest = T)
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleDResults_UniformVC", WelfareResults, envir=.GlobalEnv)
}



#################################################################
#### QV Rule D (Empirical Vote Class) - Cast Vote on Highest ####
################################################################%

# Want to choose a random class (using the empirical distribution of vote classes) and then allocate votes on highest ranked refs.

#--- determine frequency of vote classes
QV1_ClassFreq <- QV_data %>% group_by(QV1_Class) %>% summarize(Count = n()) %>% arrange(QV1_Class)
QV1_ClassFreq <- QV1_ClassFreq$Count

QV2_ClassFreq <- QV_data %>% group_by(QV2_Class) %>% summarize(Count = n()) %>% arrange(QV2_Class)
QV2_ClassFreq <- QV2_ClassFreq$Count


ruleD_EmpiricalVoteClass <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleDResults
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    # QV1
    results <- castRand(expData[row,], validRefs, ruleFreq = QV1_ClassFreq, castOnHighest = T)
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castRand(expData[row,], validRefs, ruleFreq = QV2_ClassFreq, castOnHighest = T)
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleDResults_EmpiricalVC", WelfareResults, envir=.GlobalEnv)
}


###################################
##### QV Rule D - Fifty/Fifty #####
##################################%

# First, with 50/50 chance choose the correct class vs choose a random class; Second, independent of the outcome in the first stage, with 50/50 chance allocate votes monotonically vs randomly. 


ruleD_FiftyFifty <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleDResults
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    randVector <- runif(4) # random vector used to decide what strategy to implement in the 50/50 system
    
    
    # QV1
    chosenVoteClass <- ifelse(randVector[1] < .5, chooseNormativeClass(expData[row,], validRefs), chooseRandClass(expData[row,], validRefs))
    if(randVector[2] < .5) 
    { 
      results <- castVotesMonotically(expData[row,], validRefs, chosenVoteClass) 
    } else {results <-  castVotesRandomly(expData[row,], validRefs, chosenVoteClass)}
    
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    
    # QV2
    chosenVoteClass <- ifelse(randVector[3] < .5, chooseNormativeClass(expData[row,], validRefs), chooseRandClass(expData[row,], validRefs))
    if(randVector[4] < .5) 
    { 
      results <- castVotesMonotically(expData[row,], validRefs, chosenVoteClass) 
    } else {results <-  castVotesRandomly(expData[row,], validRefs, chosenVoteClass)}
    
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleDResults_FiftyFifty", WelfareResults, envir=.GlobalEnv)
}





###########################################################
#### QV Rule E (Uniform Vote Class) - Cast Random Vote ####
##########################################################%

# Want to choose a random class (using a uniform distribution over vote classes) and then allocate votes randomly. 

ruleE_UniformVoteClass <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleEResults
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    # QV1
    results <- castRand(expData[row,], validRefs, ruleFreq = c(1,1,1,1), castOnHighest = F)
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castRand(expData[row,], validRefs, ruleFreq = c(1,1,1,1), castOnHighest = F)
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleEResults_UniformVC", WelfareResults, envir=.GlobalEnv)
}



#############################################################
#### QV Rule E (Empirical Vote Class) - Cast Random Vote ####
############################################################%

# Want to choose a random class (using the empirical distribution of vote classes) and then allocate votes randomly.  

#--- determine frequency of vote classes
QV1_ClassFreq <- QV_data %>% group_by(QV1_Class) %>% summarize(Count = n()) %>% arrange(QV1_Class)
QV1_ClassFreq <- QV1_ClassFreq$Count

QV2_ClassFreq <- QV_data %>% group_by(QV2_Class) %>% summarize(Count = n()) %>% arrange(QV2_Class)
QV2_ClassFreq <- QV2_ClassFreq$Count


ruleE_EmpiricalVoteClass <- function(expData, WelfareResults, bootstrapNum)
{
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- QV_RuleEResults
  # bootstrapNum <- 1
  
  # clear all the existing vote data
  expData <- expData %>% mutate(QV1_Class = NA, QV1_Vote_BilingualEduc = NA, QV1_Vote_Immigration = NA, QV1_Vote_PublicBonds = NA, QV1_Vote_TeacherTenure = NA,
                                QV2_Class = NA, QV2_Vote_BilingualEduc = NA, QV2_Vote_Immigration = NA, QV2_Vote_PublicBonds = NA, QV2_Vote_TeacherTenure = NA)
  
  # Implement the strategy
  for(row in 1:nrow(expData))
  {
    #--- QV1/QV2 Class and Vote Assignment
    validRefs <- refsNotAbstained(expData[row,])
    
    # QV1
    results <- castRand(expData[row,], validRefs, ruleFreq = QV1_ClassFreq, castOnHighest = F)
    expData$QV1_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV1_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
    # QV2
    results <- castRand(expData[row,], validRefs, ruleFreq = QV2_ClassFreq, castOnHighest = F)
    expData$QV2_Class[row] <- results[[1]]
    results <- detVoteVector(results[[2]])
    for(issue in c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")) {expData[[paste("QV2_Vote_",issue,sep="")]][row] <- results[[issue]][1]}
    
  }
  
  # get summaries and store them
  voteSummary <- QV_VoteSummary(expData)
  welfareSummary <- QV_WelfareSummary(expData, QV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("QV_RuleEResults_EmpiricalVC", WelfareResults, envir=.GlobalEnv)
}



###################################
###### QV Bootstrap Plot Fns ######
##################################%

# summarizes the preferences for a given proposal and puts it in a hist friendly format
# also calculates other preference metrics, e.g., num and total points in favor or against

bootPrefHist <- function(bootstrapData, Proposal)
{
  #---  Creates list of bins for histogram of preferences
  bins_width10 <- c("[-100,-90)", "[-90,-80)", "[-80,-70)", "[-70,-60)", "[-60,-50)", "[-50,-40)", 
                    "[-40,-30)", "[-30,-20)", "[-20,-10)", "[-10,0)", "0", "(0,10]", "(10,20]", 
                    "(20,30]", "(30,40]", "(40,50]", "(50,60]", "(60,70]", "(70,80]", "(80,90]", "(90,100]")
  
  bins10 <- data.frame(Bins = bins_width10, BinOrder = seq(1,21,1), stringsAsFactors = F)
  
  #--- Determines the current proposal
  
  if(Proposal == "BilingualEduc")
  {
    propTitle <- "Bilingual Education"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["EducationPref"]]
  }
  if(Proposal == "Immigration")
  {
    propTitle <- "Immigration"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["ImmigrationPref"]]
  }
  if(Proposal == "TeacherTenure")
  {
    propTitle <- "Teacher Tenure"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["TeachersPref"]]
  }
  if(Proposal == "PublicBonds")
  {
    propTitle <- "Public Bonds"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["BondsPref"]]
  }
  
  ProposalPos <- filter(bootstrapData, CurrentPref > 0) %>% mutate(CurrentPref_Bins = cut(CurrentPref, right = T, breaks = seq(0, 100, 10)))
  ProposalNeg <- filter(bootstrapData, CurrentPref < 0) %>% mutate(CurrentPref_Bins = cut(CurrentPref, right = F, breaks = seq(-100, 0, 10)))
  ProposalZero <- filter(bootstrapData, CurrentPref == 0) %>% mutate(CurrentPref_Bins = "0")
  
  Proposal <- bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, CurrentPref, CurrentPref_Bins)
  Proposal <- group_by(Proposal, CurrentPref_Bins) %>% summarise(NumSubjects = n()) %>% rename(Bins = CurrentPref_Bins)
  Proposal <- full_join(Proposal, bins10, by = "Bins")
  Proposal$Bins <- factor(Proposal$Bins, levels = Proposal$Bins[order(Proposal$BinOrder)])
  
  Proposal[is.na(Proposal)] <- 0
  Proposal <- arrange(Proposal, BinOrder) 
  
  results <- list()
  results[[1]] <- as.data.frame(t(Proposal$NumSubjects))
  results[[2]] <- sum(ProposalPos$CurrentPref)
  results[[3]] <- nrow(ProposalPos)
  results[[4]] <- abs(sum(ProposalNeg$CurrentPref))
  results[[5]] <- nrow(ProposalNeg)
  
  return(results)
}


#############################################
##### Retrieve Pop Imm Pref Data #####
############################################%

#--- Get distribution of Imm prefs from SV

# setup the binning dataframe 

bins_width10 <- c("[-100,-90)", "[-90,-80)", "[-80,-70)", "[-70,-60)", "[-60,-50)", "[-50,-40)", 
                  "[-40,-30)", "[-30,-20)", "[-20,-10)", "[-10,0)", "0", "(0,10]", "(10,20]", 
                  "(20,30]", "(30,40]", "(40,50]", "(50,60]", "(60,70]", "(70,80]", "(80,90]", "(90,100]")

bins10 <- data.frame(Bins = bins_width10, BinOrder = seq(1,21,1), stringsAsFactors = F)
bins10 <- mutate(bins10, BinColor = ifelse(Bins=="0",1,0))
binWidth <- 10

# bin and combine population data 

Pop_data <- bind_rows(select(SV_data, MID, matches("Ref")), select(QV_data, MID, matches("Ref")))
# Pop_data <- select(SV_data, MID, matches("Ref")) #***#


ProposalPos <- filter(Pop_data, ImmigrationPref > 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = T, breaks = seq(0, 100, binWidth)))
ProposalNeg <- filter(Pop_data, ImmigrationPref < 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = F, breaks = seq(-100, 0, binWidth)))
ProposalZero <- filter(Pop_data, ImmigrationPref == 0) %>% mutate(ImmigrationPref_Bins = "0")

CurrentProp <- bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, ImmigrationPref, ImmigrationPref_Bins)
CurrentProp <- group_by(CurrentProp, ImmigrationPref_Bins) %>% summarise(NumSubjects = n()) %>% rename(Bins = ImmigrationPref_Bins)
CurrentProp <- full_join(CurrentProp, bins10, by = "Bins") %>% arrange(BinOrder)
CurrentProp$Bins <- factor(CurrentProp$Bins, levels = CurrentProp$Bins[order(CurrentProp$BinOrder)])

# Create vectors with Pop immigration data
ImmPref_Pop <- CurrentProp
ImmPref_Pop[is.na(ImmPref_Pop)] <- 0
binList <- as.character(ImmPref_Pop$Bins) #***#
probDistList <- ImmPref_Pop$NumSubjects #***#

#--- Get distribution of Imm prefs from QV sample 

ProposalPos <- filter(QV_data, ImmigrationPref > 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = T, breaks = seq(0, 100, binWidth)))
ProposalNeg <- filter(QV_data, ImmigrationPref < 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = F, breaks = seq(-100, 0, binWidth)))
ProposalZero <- filter(QV_data, ImmigrationPref == 0) %>% mutate(ImmigrationPref_Bins = "0")

CurrentProp <- bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, ImmigrationPref, ImmigrationPref_Bins)
CurrentProp <- group_by(CurrentProp, ImmigrationPref_Bins) %>% summarise(NumSubjects = n()) %>% rename(Bins = ImmigrationPref_Bins)
CurrentProp <- full_join(CurrentProp, bins10, by = "Bins") %>% arrange(BinOrder)
CurrentProp$Bins <- factor(CurrentProp$Bins, levels = CurrentProp$Bins[order(CurrentProp$BinOrder)])

# Create vectors with QV immigration data and MIDs
ImmPref_QV <- CurrentProp
ImmPref_QV[is.na(ImmPref_QV)] <- 0
QV_ImmIDs <-  bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, ImmigrationPref_Bins)



###################################
###### QV Bootstrapping Loop ######
###################################

# This function will sample with replacement (boostrapping) and use the rules above to cast votes
bootstrapIterations <- 10000
set.seed(20)
sampleSize <- nrow(QV_data)

# Creates dataframes to store the Rule A-D results

QV_RuleAResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                              QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                              QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                              QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)

QV_RuleBResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                              QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                              QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                              QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)

QV_RuleBResults_classAndVote <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                                           NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                                           Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                                           QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                                           QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                                           QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                                           QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                                           AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)


QV_RuleCResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                              QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                              QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                              QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)

QV_RuleDResults_UniformVC <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                                        NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                                        Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                                        QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                                        QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                                        QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                                        QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                                        AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)

QV_RuleDResults_EmpiricalVC <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                                          NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                                          Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                                          QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                                          QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                                          QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                                          QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                                          AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)


QV_RuleDResults_FiftyFifty <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                                         NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                                         Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                                         QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                                         QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                                         QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                                         QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                                         AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)


QV_RuleEResults_UniformVC <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                                        NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                                        Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                                        QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                                        QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                                        QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                                        QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                                        AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)

QV_RuleEResults_EmpiricalVC <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                                          NQV_Result_BilingualEduc = NA, NQV_Result_Immigration = NA, NQV_Result_PublicBonds = NA, NQV_Result_TeacherTenure = NA,
                                          Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                                          QV1_Result_BilingualEduc = NA, QV1_Result_Immigration = NA, QV1_Result_PublicBonds = NA, QV1_Result_TeacherTenure = NA, 
                                          QV2_Result_BilingualEduc = NA, QV2_Result_Immigration = NA, QV2_Result_PublicBonds = NA, QV2_Result_TeacherTenure = NA,
                                          QV1_Winner_BilingualEduc = NA, QV1_Winner_Immigration = NA, QV1_Winner_PublicBonds = NA, QV1_Winner_TeacherTenure = NA, 
                                          QV2_Winner_BilingualEduc = NA, QV2_Winner_Immigration = NA, QV2_Winner_PublicBonds = NA, QV2_Winner_TeacherTenure = NA,
                                          AggregateWelfare_noQV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_QV1 = NA, AggregateWelfare_QV2 = NA)


# Creates dataframes to store ALL preferences for all proposals and bootstraps.
QV_EducBSPref <- data.frame()
QV_BondBSPref <- data.frame()
QV_TeachBSPref <- data.frame()
QV_ImmBSPref <- data.frame()

QV_PrefSummary <- data.frame(Round = seq(1, bootstrapIterations, 1),
                             BilingualEduc_PointsInFavor = NA, BilingualEduc_SubjectsInFavor = NA, 
                             BilingualEduc_PointsAgainst = NA, BilingualEduc_SubjectsAgainst = NA, 
                             Immigration_PointsInFavor = NA, Immigration_SubjectsInFavor = NA, 
                             Immigration_PointsAgainst = NA, Immigration_SubjectsAgainst = NA,
                             PublicBonds_PointsInFavor = NA, PublicBonds_SubjectsInFavor = NA,
                             PublicBonds_PointsAgainst = NA, PublicBonds_SubjectsAgainst = NA,
                             TeacherTenure_PointsInFavor = NA, TeacherTenure_SubjectsInFavor = NA, 
                             TeacherTenure_PointsAgainst = NA, TeacherTenure_SubjectsAgainst = NA)

# Creates dataframes to store the realized individual utilities in each bootstrap
QV_IndivUtility_NQV_A <- data.frame()
QV_IndivUtility_QV1_A <- data.frame()
QV_IndivUtility_QV2_A <- data.frame()

QV_IndivUtility_NQV_B <- data.frame()
QV_IndivUtility_QV1_B <- data.frame()
QV_IndivUtility_QV2_B <- data.frame()

QV_IndivUtility_NQV_B_classAndVote <- data.frame()
QV_IndivUtility_QV1_B_classAndVote <- data.frame()
QV_IndivUtility_QV2_B_classAndVote <- data.frame()

QV_IndivUtility_NQV_C <- data.frame()
QV_IndivUtility_QV1_C <- data.frame()
QV_IndivUtility_QV2_C <- data.frame()

QV_IndivUtility_NQV_D_Emp <- data.frame()
QV_IndivUtility_QV1_D_Emp <- data.frame()
QV_IndivUtility_QV2_D_Emp <- data.frame()

QV_IndivUtility_NQV_D_Unif <- data.frame()
QV_IndivUtility_QV1_D_Unif <- data.frame()
QV_IndivUtility_QV2_D_Unif <- data.frame()

QV_IndivUtility_NQV_D_5050 <- data.frame()
QV_IndivUtility_QV1_D_5050 <- data.frame()
QV_IndivUtility_QV2_D_5050 <- data.frame()

QV_IndivUtility_NQV_E_Emp <- data.frame()
QV_IndivUtility_QV1_E_Emp <- data.frame()
QV_IndivUtility_QV2_E_Emp <- data.frame()

QV_IndivUtility_NQV_E_Unif <- data.frame()
QV_IndivUtility_QV1_E_Unif <- data.frame()
QV_IndivUtility_QV2_E_Unif <- data.frame()

# function to datermine indiv util vectors under diff voting schemes

calcRealizedUtil <- function(BS_data, QV_summary)
{

  QV_temp <- BS_data %>% select(-matches("Abstain")) %>% select(matches("Ref")) %>%
    mutate(Util_NQV_BilingualEduc = ifelse(Ref1_BilingualEduc == QV_summary$NoQV_Result[QV_summary$Referenda == "BilingualEduc"],EducationPref,0),
           Util_QV1_BilingualEduc = ifelse(Ref1_BilingualEduc == QV_summary$QV1_Result[QV_summary$Referenda == "BilingualEduc"],EducationPref,0),
           Util_QV2_BilingualEduc = ifelse(Ref1_BilingualEduc == QV_summary$QV2_Result[QV_summary$Referenda == "BilingualEduc"],EducationPref,0),
           Util_NQV_TeacherTenure = ifelse(Ref2_TeacherTenure == QV_summary$NoQV_Result[QV_summary$Referenda == "TeacherTenure"],TeachersPref,0),
           Util_QV1_TeacherTenure = ifelse(Ref2_TeacherTenure == QV_summary$QV1_Result[QV_summary$Referenda == "TeacherTenure"],TeachersPref,0),
           Util_QV2_TeacherTenure = ifelse(Ref2_TeacherTenure == QV_summary$QV2_Result[QV_summary$Referenda == "TeacherTenure"],TeachersPref,0),
           Util_NQV_PublicBonds = ifelse(Ref3_PublicBonds == QV_summary$NoQV_Result[QV_summary$Referenda == "PublicBonds"],BondsPref,0),
           Util_QV1_PublicBonds = ifelse(Ref3_PublicBonds == QV_summary$QV1_Result[QV_summary$Referenda == "PublicBonds"],BondsPref,0),
           Util_QV2_PublicBonds = ifelse(Ref3_PublicBonds == QV_summary$QV2_Result[QV_summary$Referenda == "PublicBonds"],BondsPref,0),
           Util_NQV_Immigration = ifelse(Ref4_Immigration == QV_summary$NoQV_Result[QV_summary$Referenda == "Immigration"],ImmigrationPref,0),
           Util_QV1_Immigration = ifelse(Ref4_Immigration == QV_summary$QV1_Result[QV_summary$Referenda == "Immigration"],ImmigrationPref,0),
           Util_QV2_Immigration = ifelse(Ref4_Immigration == QV_summary$QV2_Result[QV_summary$Referenda == "Immigration"],ImmigrationPref,0)) %>%
    mutate(Util_NQV = Util_NQV_BilingualEduc + Util_NQV_TeacherTenure + Util_NQV_PublicBonds + Util_NQV_Immigration,
           Util_QV1 = Util_QV1_BilingualEduc + Util_QV1_TeacherTenure + Util_QV1_PublicBonds + Util_QV1_Immigration, 
           Util_QV2 = Util_QV2_BilingualEduc + Util_QV2_TeacherTenure + Util_QV2_PublicBonds + Util_QV2_Immigration)
  
  temp_NQV <- as.data.frame(t(QV_temp$Util_NQV))
  temp_QV1 <- as.data.frame(t(QV_temp$Util_QV1))
  temp_QV2 <- as.data.frame(t(QV_temp$Util_QV2))
  
  return(list(temp_NQV,temp_QV1,temp_QV2))
} 




# where each simulated dataframe will be stored

simulatedSamples <- list()


# lists to save summary of each simulation data analysis
list_Summary_Pre <- data.frame()
list_Summary_Mid <- data.frame()
list_Summary_Final <- data.frame()


# This loop bootstraps a new sample for QV of the same size using sampling w/ replacement
for(bsCount in 1:bootstrapIterations) 
{
  print(bsCount)
  
  BS_data <- data.frame()
  
  # Use QV_Imm data to create a random sample and then draw from QV_ImmIDs
  randomBinSample <- sample(binList, nrow(QV_data), replace = T, prob = probDistList)
  randomBinSample <- as.data.frame(table(randomBinSample))
  randomBinSample <- rename(randomBinSample, Bin = randomBinSample) %>% mutate(Bin = as.character(Bin))
  
  for(row in 1:nrow(randomBinSample))
  {
    # Subset the list of names by bin
    QV_temp <- filter(QV_ImmIDs, ImmigrationPref_Bins == randomBinSample$Bin[row])
    
    # Take a random sample and store only the MIDs
    randRowNum <- sample(nrow(QV_temp), randomBinSample$Freq[row], replace = T)
    QV_temp <- QV_temp[randRowNum,]
    QV_temp <- select(QV_temp, MID)
    
    # Merge random MID sample with final MID list
    BS_data <- bind_rows(BS_data, QV_temp)
  }
  
  # Use final bootstrapped list to create a new QV sample
  selectedRows <- match(BS_data$MID, QV_data$MID)
  BS_data <- QV_data[selectedRows,]
  
  
  # Number the rows to avoid conflicts
  row.names(BS_data) <- seq(1,nrow(BS_data),1)
  
  #save an unedited copy for the preference summary
  BS_data_temp <- BS_data
  simulatedSamples[[bsCount]] <- BS_data
  
  
  # Mutate the data back to all positive & Run the bootstrapped data under the four rules
  BS_data <- mutate(BS_data, ImmigrationPref = abs(ImmigrationPref), BondsPref = abs(BondsPref),
                    EducationPref = abs(EducationPref), TeachersPref = abs(TeachersPref))
  
  # Testing
  # expData <- BS_data
  # WelfareResults <- QV_RuleAResults
  # bootstrapNum <- bsCount
  
  # Run the bootstrapped data under the rules
  
  temp_Pre_All <- data.frame()
  temp_Mid_All <- data.frame()
  temp_Final_All <- data.frame()
  
  ruleA(expData = BS_data, WelfareResults = QV_RuleAResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "A", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "A", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "A", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_A <- bind_rows(QV_IndivUtility_NQV_A, utils[[1]])
  QV_IndivUtility_QV1_A <- bind_rows(QV_IndivUtility_QV1_A, utils[[2]]) 
  QV_IndivUtility_QV2_A <- bind_rows(QV_IndivUtility_QV2_A, utils[[3]])
  
  ruleB(expData = BS_data, WelfareResults = QV_RuleBResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "B", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "B", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "B", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_B <- bind_rows(QV_IndivUtility_NQV_B, utils[[1]])
  QV_IndivUtility_QV1_B <- bind_rows(QV_IndivUtility_QV1_B, utils[[2]]) 
  QV_IndivUtility_QV2_B <- bind_rows(QV_IndivUtility_QV2_B, utils[[3]])
  
  ruleB_classAndVote(expData = BS_data, WelfareResults = QV_RuleBResults_classAndVote, bootstrapNum = bsCount)

  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "B - CNV", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "B - CNV", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "B - CNV", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_B_classAndVote <- bind_rows(QV_IndivUtility_NQV_B_classAndVote, utils[[1]])
  QV_IndivUtility_QV1_B_classAndVote <- bind_rows(QV_IndivUtility_QV1_B_classAndVote, utils[[2]])
  QV_IndivUtility_QV2_B_classAndVote <- bind_rows(QV_IndivUtility_QV2_B_classAndVote, utils[[3]])

  
  ruleC(expData = BS_data, WelfareResults = QV_RuleCResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "C", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "C", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "C", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_C <- bind_rows(QV_IndivUtility_NQV_C, utils[[1]])
  QV_IndivUtility_QV1_C <- bind_rows(QV_IndivUtility_QV1_C, utils[[2]]) 
  QV_IndivUtility_QV2_C <- bind_rows(QV_IndivUtility_QV2_C, utils[[3]])
  
  ruleD_EmpiricalVoteClass(expData = BS_data, WelfareResults = QV_RuleDResults_EmpiricalVC, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "D - Emp", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "D - Emp", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "D - Emp", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_D_Emp <- bind_rows(QV_IndivUtility_NQV_D_Emp, utils[[1]])
  QV_IndivUtility_QV1_D_Emp <- bind_rows(QV_IndivUtility_QV1_D_Emp, utils[[2]]) 
  QV_IndivUtility_QV2_D_Emp <- bind_rows(QV_IndivUtility_QV2_D_Emp, utils[[3]])
  
  ruleD_UniformVoteClass(expData = BS_data, WelfareResults = QV_RuleDResults_UniformVC, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "D - Unif", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "D - Unif", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "D - Unif", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_D_Unif <- bind_rows(QV_IndivUtility_NQV_D_Unif, utils[[1]])
  QV_IndivUtility_QV1_D_Unif <- bind_rows(QV_IndivUtility_QV1_D_Unif, utils[[2]]) 
  QV_IndivUtility_QV2_D_Unif <- bind_rows(QV_IndivUtility_QV2_D_Unif, utils[[3]])
  
  ruleD_FiftyFifty(expData = BS_data, WelfareResults = QV_RuleDResults_FiftyFifty, bootstrapNum = bsCount)

  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "D - 50/50", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "D - 50/50", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "D - 50/50", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_D_5050 <- bind_rows(QV_IndivUtility_NQV_D_5050, utils[[1]])
  QV_IndivUtility_QV1_D_5050 <- bind_rows(QV_IndivUtility_QV1_D_5050, utils[[2]])
  QV_IndivUtility_QV2_D_5050 <- bind_rows(QV_IndivUtility_QV2_D_5050, utils[[3]])
  
  ruleE_EmpiricalVoteClass(expData = BS_data, WelfareResults = QV_RuleEResults_EmpiricalVC, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "E - Emp", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "E - Emp", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "E - Emp", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_E_Emp <- bind_rows(QV_IndivUtility_NQV_E_Emp, utils[[1]])
  QV_IndivUtility_QV1_E_Emp <- bind_rows(QV_IndivUtility_QV1_E_Emp, utils[[2]]) 
  QV_IndivUtility_QV2_E_Emp <- bind_rows(QV_IndivUtility_QV2_E_Emp, utils[[3]])
  
  ruleE_UniformVoteClass(expData = BS_data, WelfareResults = QV_RuleEResults_UniformVC, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "E - Unif", Testing_QV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "E - Unif", Testing_QV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "E - Unif", QV_summary))
  utils <- calcRealizedUtil(BS_data, QV_summary)
  QV_IndivUtility_NQV_E_Unif <- bind_rows(QV_IndivUtility_NQV_E_Unif, utils[[1]])
  QV_IndivUtility_QV1_E_Unif <- bind_rows(QV_IndivUtility_QV1_E_Unif, utils[[2]]) 
  QV_IndivUtility_QV2_E_Unif <- bind_rows(QV_IndivUtility_QV2_E_Unif, utils[[3]])
  
  # combine and save data 
  
  list_Summary_Pre <- rbind(list_Summary_Pre, temp_Pre_All)
  list_Summary_Mid <- rbind(list_Summary_Mid, temp_Mid_All)
  list_Summary_Final <- rbind(list_Summary_Final, temp_Final_All)
  
  
  
  #--- Preference Histogram and other Preference Summaries 
  
  # BilingualEduc
  PrefHistResults <- bootPrefHist(BS_data_temp, "BilingualEduc")
  QV_EducBSPref <- bind_rows(QV_EducBSPref, PrefHistResults[[1]])
  
  QV_PrefSummary$BilingualEduc_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  QV_PrefSummary$BilingualEduc_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  QV_PrefSummary$BilingualEduc_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  QV_PrefSummary$BilingualEduc_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
  
  # Immigration 
  PrefHistResults <- bootPrefHist(BS_data_temp, "Immigration")
  QV_ImmBSPref <- bind_rows(QV_ImmBSPref, PrefHistResults[[1]])
  
  QV_PrefSummary$Immigration_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  QV_PrefSummary$Immigration_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  QV_PrefSummary$Immigration_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  QV_PrefSummary$Immigration_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
  # TeacherTenure 
  PrefHistResults <- bootPrefHist(BS_data_temp, "TeacherTenure")
  QV_TeachBSPref <- bind_rows(QV_TeachBSPref, PrefHistResults[[1]])
  
  QV_PrefSummary$TeacherTenure_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  QV_PrefSummary$TeacherTenure_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  QV_PrefSummary$TeacherTenure_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  QV_PrefSummary$TeacherTenure_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
  # PublicBonds 
  PrefHistResults <- bootPrefHist(BS_data_temp, "PublicBonds")
  QV_BondBSPref <- bind_rows(QV_BondBSPref, PrefHistResults[[1]])
  
  QV_PrefSummary$PublicBonds_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  QV_PrefSummary$PublicBonds_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  QV_PrefSummary$PublicBonds_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  QV_PrefSummary$PublicBonds_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
}



# Rename the column names of the dataframes with the preferences for each bootstrap and each referenda
names(QV_EducBSPref) <- gsub("V","EducPref",names(QV_EducBSPref))
names(QV_ImmBSPref) <- gsub("V","ImmPref",names(QV_ImmBSPref))
names(QV_BondBSPref) <- gsub("V","BondPref",names(QV_BondBSPref))
names(QV_TeachBSPref) <- gsub("V","TeachPref",names(QV_TeachBSPref))


# Create a column to store the Round number for the list of individual utilities for each bootstrap
QV_IndivUtility_NQV_A <- tibble::rownames_to_column(QV_IndivUtility_NQV_A, "Round")
QV_IndivUtility_NQV_B <- tibble::rownames_to_column(QV_IndivUtility_NQV_B, "Round")
QV_IndivUtility_NQV_B_classAndVote <- tibble::rownames_to_column(QV_IndivUtility_NQV_B_classAndVote, "Round")
QV_IndivUtility_NQV_C <- tibble::rownames_to_column(QV_IndivUtility_NQV_C, "Round")
QV_IndivUtility_NQV_D_5050 <- tibble::rownames_to_column(QV_IndivUtility_NQV_D_5050, "Round")
QV_IndivUtility_NQV_D_Unif <- tibble::rownames_to_column(QV_IndivUtility_NQV_D_Unif, "Round")
QV_IndivUtility_NQV_D_Emp <- tibble::rownames_to_column(QV_IndivUtility_NQV_D_Emp, "Round")
QV_IndivUtility_NQV_E_Unif <- tibble::rownames_to_column(QV_IndivUtility_NQV_E_Unif, "Round")
QV_IndivUtility_NQV_E_Emp <- tibble::rownames_to_column(QV_IndivUtility_NQV_E_Emp, "Round")
QV_IndivUtility_QV1_A <- tibble::rownames_to_column(QV_IndivUtility_QV1_A, "Round")
QV_IndivUtility_QV2_A <- tibble::rownames_to_column(QV_IndivUtility_QV2_A, "Round")
QV_IndivUtility_QV1_B <- tibble::rownames_to_column(QV_IndivUtility_QV1_B, "Round")
QV_IndivUtility_QV2_B <- tibble::rownames_to_column(QV_IndivUtility_QV2_B, "Round")
QV_IndivUtility_QV1_B_classAndVote <- tibble::rownames_to_column(QV_IndivUtility_QV1_B_classAndVote, "Round")
QV_IndivUtility_QV2_B_classAndVote <- tibble::rownames_to_column(QV_IndivUtility_QV2_B_classAndVote, "Round")
QV_IndivUtility_QV1_C <- tibble::rownames_to_column(QV_IndivUtility_QV1_C, "Round")
QV_IndivUtility_QV2_C <- tibble::rownames_to_column(QV_IndivUtility_QV2_C, "Round")
QV_IndivUtility_QV1_D_5050 <- tibble::rownames_to_column(QV_IndivUtility_QV1_D_5050, "Round")
QV_IndivUtility_QV2_D_5050 <- tibble::rownames_to_column(QV_IndivUtility_QV2_D_5050, "Round")
QV_IndivUtility_QV1_D_Unif <- tibble::rownames_to_column(QV_IndivUtility_QV1_D_Unif, "Round")
QV_IndivUtility_QV2_D_Unif <- tibble::rownames_to_column(QV_IndivUtility_QV2_D_Unif, "Round")
QV_IndivUtility_QV1_D_Emp <- tibble::rownames_to_column(QV_IndivUtility_QV1_D_Emp, "Round")
QV_IndivUtility_QV2_D_Emp <- tibble::rownames_to_column(QV_IndivUtility_QV2_D_Emp, "Round")
QV_IndivUtility_QV1_E_Unif <- tibble::rownames_to_column(QV_IndivUtility_QV1_E_Unif, "Round")
QV_IndivUtility_QV2_E_Unif <- tibble::rownames_to_column(QV_IndivUtility_QV2_E_Unif, "Round")
QV_IndivUtility_QV1_E_Emp <- tibble::rownames_to_column(QV_IndivUtility_QV1_E_Emp, "Round")
QV_IndivUtility_QV2_E_Emp <- tibble::rownames_to_column(QV_IndivUtility_QV2_E_Emp, "Round")



# Save and export the dataframes with the welfare results and bootsrapped preferences 
write.csv(QV_EducBSPref, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_EducBSPref.csv", row.names = F)
write.csv(QV_ImmBSPref, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_ImmBSPref.csv", row.names = F)
write.csv(QV_BondBSPref, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_BondBSPref.csv", row.names = F)
write.csv(QV_TeachBSPref, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_TeachBSPref.csv", row.names = F)

write.csv(QV_PrefSummary, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_PrefPointSummary.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_A, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_A.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_B, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_B.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_B_classAndVote, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_B_CNV.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_C, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_C.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_D_Emp, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_D_Emp.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_D_Unif, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_D_Unif.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_D_5050, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_D_5050.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_E_Emp, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_E_Emp.csv", row.names = F)
write.csv(QV_IndivUtility_NQV_E_Unif, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_E_Unif.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_A, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_A.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_A, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_A.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_B, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_B.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_B, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_B.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_B_classAndVote, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_B_CNV.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_B_classAndVote, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_B_CNV.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_C, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_C.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_C, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_C.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_D_Emp, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_D_Emp.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_D_Emp, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_D_Emp.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_D_Unif, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_D_Unif.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_D_Unif, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_D_Unif.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_D_5050, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_D_5050.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_D_5050, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_D_5050.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_E_Emp, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_E_Emp.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_E_Emp, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_E_Emp.csv", row.names = F)
write.csv(QV_IndivUtility_QV1_E_Unif, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_E_Unif.csv", row.names = F)
write.csv(QV_IndivUtility_QV2_E_Unif, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_E_Unif.csv", row.names = F)

write.csv(list_Summary_Pre, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_Summary_Pre.csv", row.names = F)
write.csv(list_Summary_Mid, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_Summary_Mid.csv", row.names = F)
write.csv(list_Summary_Final, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_Summary_Final.csv", row.names = F)

write.csv(QV_RuleAResults, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleAResults.csv", row.names = F)
write.csv(QV_RuleBResults, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleBResults.csv", row.names = F)
write.csv(QV_RuleBResults_classAndVote, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleBResults_CNV.csv", row.names = F)
write.csv(QV_RuleCResults, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleCResults.csv", row.names = F)
write.csv(QV_RuleDResults_EmpiricalVC, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleDResults_EmpiricalVC.csv", row.names = F)
write.csv(QV_RuleDResults_UniformVC, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleDResults_UniformVC.csv", row.names = F)
write.csv(QV_RuleDResults_FiftyFifty, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleDResults_FiftyFifty.csv", row.names = F)
write.csv(QV_RuleEResults_EmpiricalVC, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleEResults_EmpiricalVC.csv", row.names = F)
write.csv(QV_RuleEResults_UniformVC, "./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleEResults_UniformVC.csv", row.names = F)

saveRDS(simulatedSamples, file="./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/simulatedSamples.RData")